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Abstract 


In this paper we consider fourth order difference approximations of initial-boundary value 
problems for hyperbolic partial differential equations. We use the method of lines ap¬ 
proach with both explicit and compact implicit difference operators in space. The explicit 
operator satisfies an energy estimate leading to strict stability. For the implicit operator 
we develop boundary conditions and give a complete proof of strong stability using the 
Laplace transform technique. 

We also present numerical experiments for the linear advection equation and Burgers’ 
equation with discontinuities in the solution or in its derivative. The first equation is used 
for modeling contact discontinuities in fluid dynamics, the second one for modeling shocks 
and rarefaction waves. The time discretization is done with a third order Runge-Kutta 

TVD method. For solutions with discontinuities in the solution itself we add a filter based 
on second order viscosity. 

In case of the non-linear Burgers’ equation we use a flux splitting technique that results 
in an energy estimate tor certain difference approximations, in which case also an entropy 
condition is fulfilled. In particular we shall demonstrate that the unsplit conservative 
form produces a non-physical shock instead of the physically correct rarefaction wave. 

In the numerical experiments we compare our fourth order methods with a standard 
second order one and with a third order TVD-method. The results show that the fourth 
order methods are the only ones that give good results for all the considered test problems. 



1 Introduction 

It is well known that high-order accurate difference operators are more efficient than low- 
order ones for hyperbolic problems with smooth solutions, except for very low accuracy 
requirements in the solution. The theoretical basis for this conclusion is found in [7] 
and [17]. Nevertheless, in practice most calculations are done with first or second order 
approximations. One of the reasons for this is the extra difficulty that arises near the 
boundaries. It is always possible to derive non-symmetric operators near the boundaries 
that have sufficient formal accuracy, but it is more difficult when requiring that the method 
also be stable. In [12] and [16] high-order methods for initial-boundary value problems 
are constructed based on the work by Kreiss and Scherer [8] and [9]. The approximations 
satisfy an energy estimate that guarantees strict stability. For integrations over long time 
intervals this is an especially important property. 

Stability an. ysis based on Laplace transform technique leads to strong stability if 
the Kreiss condition is satisfied as shown in [5]. In this book there is also a complete 
analysis of a semi-discrete fourth order approximation based on the standard five-point 
scheme, and the Kreiss condition is shown to be satisfied. Strict stability, however, is not 
an automatic consequence of this theory. 

In sec. 2 we give a brief review of the currently available results for fourth order 
accurate operators. 

Compact difference operators (Pade type) for the space part of PDEs have been con¬ 
sidered, for example, in [17], [14] and [10]. These methods are based on an approximation 

—► P _1 Q, where P and Q are non-diagonal difference operators. In this way the error 
constant can be substantially reduced, and the extra work required for solving the banded 
systems in each time step may well pay off. 

In [1] a boundary procedure is developed for the fourth order case, where P and Q are 
tridiagonal, and it is verified that the Kreiss condition is satisfied. However, the step from 
the Kreiss condition to stability is not carried out. No such general theory is currently 
available; in [5] only the explicit case P = 1 is treated. In sec. 3 we present the full 
theory for the implicit fourth order approximation by generalizing the Laplace transform 
technique. We construct boundary conditions such that the resulting approximation is 
strongly stable and gives a fourth order error estimate. 

For problems with non-smooth solutions, the error estimates based on the truncation 
error breaks down. Still the fourth order methods may well also be competitive with lower 
order ones in this case. This is demonstrated in sec. 4, where we present a number of 
numerical experiments. We use the linear advection equation and Burgers’ equation with 
discontinuities in the solution or in its derivative. The first equation is a good model for 
contact discontinuities in fluid dynamics; the second one is used for modeling shocks and 
rarefaction waves. The time discretization is done with a third order Runge-Kutta TVD 
method. 
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For solutions with discontinuities in the derivatives, for example rarefaction waves, no 
extra viscosity terms are necessary. However, for discontinuities in the solution itself, we 
expect oscillations in the numerical solution. Therefore, we add a filter based on second 
order viscosity. This takes the formal accuracy down to first order, but by using a switch 
as coefficient for the viscosity term, this loss of accuracy is limited to the immediate area 
around the discontinuity. 

In case of the non-linear Burgers' equation we use a flux splitting technique that results 
in an energy estimate for certain difference approximations, in which case also an entropy 
condition is fulfilled as described in [13J. In particular we shall demonstrate that the 
unsplit conservative form produces a non-physical shock instead of a rarefaction wave. 


2 Explicit Difference Operators 


It is common to use the energy method in order to establish well-posedness of initial¬ 
boundary value problems (IBVP) for hyperbolic PDEs. Consider the model problem 


ti* 4* u x = 0, 

“(M) = <7(0. 
u(i,0) = /(*), 


0 < x < oo, / > 0 


0 ) 


where the initial data f(x) is assumed to have compact support. We consider the quarter 
space problem for convenience; domains with two boundaries are handled analogously, 
cf. sec. 4. The standard scalar product and the corresponding norm are defined by 

roo 

(u,v)= u(x)v(x)dx, j|u|| 2 = (u,u). 

JO 

To arrive at an a priori estimate for eq. (1) we use the following tools. 

(i) Integration by parts (assuming compact support): 


d 

dt 


— Hull 2 = 2 (u,u t ) = -2(t/,u r ) = u(0,0 2 • 


(ii) Boundary conditions: 


Hence, 


«(0,0 = g{t)- 


Jt M 2 = g{t)\ 

which after integration with respect to t yields an energy estimate 

IW-. Oil’ = ll/ll 1 +/‘l»(r)|Vr. 

Jo 
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For the outflow problem 


tit — u x = 0, 0<x<oo. t>0 

u(x,0 )*/(*), (2) 

no boundary condition is needed to obtain an energy estimate; in fact, one can estimate 
the solution at the boundary x = 0: 

/'|u(0,r)| ! dr = ll/H*. 

JO 

It is also possible to derive an energy estimate for the nonlinear conservation law 

u t + F x = 0, 0 < x < oo , t > 0, 

«(M)®$(0 ifF'(u)>0, (3) 

«(x,0) = /(x), 

provided F(u) satisfies a certain structural hypothesis. The key to obtaining an energy 
estimate lies in splitting the flux derivative F x into two parts, 

Fr = (F - G) x + G x = (F - G) t + G'u z , 

where G — G(u) satisfies Euler’s inhomogeneous differential equation 

G'u = -G + F <=» G = - T F(v)dv. 

u Jo 

Hence, F x can be written as 

F x - (G'u) x + G'u x , 

which will be referred to as the canonical splitting of F x . The solution of eq. (3) then 
satisfies 

d 

^||u|| 2 = —2(u, F x ) = -2(u,(G , u) x )-2(u,G'u r ) = 2uC? / u(0,<), 

where 

uG'u— f F'{v)vdv. (4) 

Jo 

Thus, in order to obtain an energy estimate we must confine ourselves to flux functions 
F such that the sign of F'(u) determines that of (4). This is true if, for instance, 

sgn(u) = sgn(F'(u)) or sgn(F'(u)) > (<) 0. 

The former condition is true for Burgers’ equation, whereas the latter holds for all linear, 
constant coefficient equations. We thus have an example of the previously mentioned 
structural hypothesis. The canonical splitting and the structural hypothesis can be gen¬ 
eralized to symmetrizabie systems and several space dimensions [13]. Hence, if we are to 
obtain an energy estimate for a nonlinear conservation law, the list of tools is augmented 
by 
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(iii) Canonical splitting of F x . 

(iv) Structural hypothesis on F. 


The analysis of the semi-discrete case can be carried out in much the same way as in 
the continuous case; integration by parts is replaced by summation by parts. The main 
difficulty lies in the treatment of the analytic boundary conditions. 

The discrete L 2 scalar product and norm are defined as 


(tl, v)o,oo — •> ll^llo.oo (tt,u)o,oo- 

;=° 

To make the notation less cumbersome we shall use the conventions (u,v) = (u,r)o,oo and 
|| U || = ||u||o,oo- The difference operator D is defined by 


m 


^ , J — 0, 1, . . . y 


k=l 


where D is a local operator, i. e., |/, - j\ < l, \m } - j\ < m for some constants /, m; h is 
the (uniform) mesh size. For certain operators one can find a local, symmetric positive 
definite operator (SPD) H [8, 9, 3, 2], such that 


(u, HDv) = -u Q v 0 - ( DxiyHv) 



in complete analogy with the analytic case. As usual we have assumed compact support. 
It can be shown [8] that it is impossible to choose H = / for consistent approximations 
D. An example of a fourth order accurate operator with third order boundary closure 
satisfying eq. (5) is provided in the Appendix. 

The treatment of the analytic boundary conditions can be done in various ways. One 
possibility is to represent the analytic boundary conditions as a projection operator T. 
Eq.(l) is then discretized as 

^ + (TDu), = «<-T) d f t ),, J —0,1 . (6) 

»,( 0 ) = !, 

where 

(Tu) 0 = 0 , {Tu) } = Uj, j = 1 , 2 ,... , 

and g — (g(t) x .. .) r ; <?(*) is the analytic boundary data, and x is a generic component. 
The actual value is of no importance. If uo(0) — fo — <?(0) it follows that uo(f) = g(t) or, 
equivalently, (I - T)(u -j) = 0. Hence, the analytic boundary condition is fulfilled for 
all time. Assuming that / and g satisfy certain compatibility conditions one can prove 
the following estimate (12) 

Hu(t)) = (/,///)+ /V(r)dr. 

Jo 
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Since H is a bounded, symmetric positive definite operator the above equation yields 

IMP < const. (ll/H 2 + £9 2 {T)dr\ . 

More generally, given a norm H, any linear boundary condition can be represented 

“ %~i OPer t, r T SUCh that HT = TTH [12 J- This P r °P ert y together with 
eq. (5) makes it possible to prove strict stability for arbitrarily accurate semi-discrete 

approximations of hyperbolic systems in one space dimension. By strict stability we 

mean that the growth rate of the analytic and the semi-discrete solutions is identical 

Confining ourselves to diagonal norms H it can be shown that operators D satisfying 

eq. (5) will result ,n strictly stable semi-discrete approximations of hyperbolic systems 

in several space dimensions. The stability results are valid for curvilinear domains with 

hi > nh S rr h d ff OUndar,eS ’ Cf ‘ [12) f ° r 4 Complete Presentation. For explicit examples of 

high-order difference operators corresponding to diagonal norms we refer to fl 1 2 161 If 

we relax eq. (5) to 1 ’ ’ 11 


(u,HDv) - B{u b ,v b ) - {Du,Hv), u b = (« 0 ... u q ) T , u 6 = (u 0 


v i) 


(7) 


for some function B and some constant q, it is in general no longer possible to prove stricl 
stability. However, ,t may still be possible to prove stability using Laplace transform 
techniques [5]; at outflow boundaries one uses extrapolation, and at inflow boundaries 
the differential equation is used to impose proper analytic boundary conditions, cf. sec. 3 
Yet another technique for enforcing the analytic boundary conditions is used [3], where a 

Tea %\ U T° n T_! ta / ne r ApPr ° Xirnati0n Term ) is a dded to the right hand side 
q. (6) after setting T - /. The penalty function is constructed such that the solution 

of the semi-discrete scheme will satisfy the analytic boundary conditions to some order 

of accuracy. It can be shown that the resulting semi-discrete scheme is strictly stable 

or one- imensional constant-coefficient hyperbolic systems. Finally we mention that the 

projection technique outlined above carries over to the nonlinear case if the semi-discrete 

equation is based on the canonical splitting of the flux derivative, and if D satisfies eq. (5) 

for some diagonal norm. This analysis will be carried out in a forthcoming paper. 

3 Implicit Difference Operators 


In this section we shall construct boundary conditions for the standard fourth orde 
implicit approximation and prove stability. It has been shown in [4] that it is impossibl 
to enforce eq. (5) for sufficiently accurate boundary approximations as long as the matri 
H in the norm is non-diagonal only in a neighborhood of the boundary. Therefore, we shal 
use the Laplace transform technique to prove stability. Note that for the two-boundar 
problem, however, the semi-discrete solution may grow exponentially in time even if th 
analytic solution is bounded in time This is in agreement with the discussion in sec 2 
since one cannot in general prove st stability using the Laplace transform technique 
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The step from the Kreiss condition to the stability estimate is not covered by existing 
theory. We shall use the same type of technique as used for explicit approximations in 
[5], but it will be modified so as to apply to implicit operators. 

We first consider the outflow problem 

u t = u x , 0 < z < oo,0 < t 

u(z,0) = f{x). 

Let vj be the approximation of u x (xj,t). The standard fourth order implicit approxima¬ 
tion used at inner points is 


+ Ay, + t>j+i) = ^-(u J+1 ~ wj-i) • j = 1,2,... . 

Since there is no boundary condition for u at x = 0, we use a one-sided approximation at 
j = 0. A Taylor expansion shows that 

Vo 4- 2v\ = —(—5«o + 4tii + Uj) 


has a truncation error of order h 3 (for a systematic derivation of high-order approxima¬ 
tions, see [10]). 

Let the operators P, Q be defined by 


t 


(Pu), = { 


\ 


— (u 0 + 2ui), 

+ 4uj +u J+1 ), 



( 1 


(Qu), = { 


24 h 
1 


(—5tio + 4uj + tij), 


y 2 h 


K+i 




where the boundary approximation has been normalized so that P is symmetric. For 
general problems, we solve for the approximation v of u z from 


(Pv), = (Qu)j, j 



and substitute v into the general approximation of the differential equation. For our 
model problem, the approximation can be written as 


(p 7t ) ’ = (Qu) '’ J=0 ’ 1 ’ 
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We need to know that our approximation is solvable. Furthermore, the operator P is 
going to be used to define a norm. We have 

Lemma 3.1 P is a symmetric positive definite operator in Co. 


Proof: The matrix representation of P shows that it is symmetric. We have 

— (u, Pu) — (|uo| 2 4- 2uoUi) + (2ujUo 4- 8|ui| 2 + 2 uiu 2 ) 
h 

oo 

+(2u 2 u, 4- 8|u 2 | 2 4- 2u 2 u 3 ) + ... > |uo| 2 - 4|u 0 «i | 4- 6|«i I 2 + 4 £ NI* 

J=2 

. oo J 

> |u 0 | 2 - -Kl 2 — 5|ui| 2 4- 6|ui| 2 4- 4 ^ |«j| 2 > ’ 

O jrs2 

which proves the lemma. 

For the purpose of deriving stability and error estimates, it is convenient to rewrite 
the approximation with the boundary scheme singled out. With inhomogeneous terms in 
the boundary approximation, (10) becomes 

(P~h = J = '- 2 . 

du (11) 

(P-fi)o = (<M)4-S, 

u j(H) = fi • 


We prove 

Lemma 3.2 Consider (It) with / = 0. The solution satisfies the estimate 
||u(0!! 2 + J |uj(r)| 2 dr < const. \hg(r)\ 2 dr , J =0,1,.... 

Proof: The Laplace transformed approximation is 

s(Pu), = (Qu)>, j = 1»2,... , 

s{Pu) 0 = (<?u)o 4 g . 

INI < 

with the characteristic equation at inner points 

i(x 2 4 4 k 4-1) = 3(/c 2 - 1), s = sh. 


( 12 ) 

(13) 

(14) 

(15) 
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This equation has exactly one root k x with |«i| < 1 for Re(s) > 0. In order to prove 
this, we first note that there is no root k with |k| = 1 for Re(5) > 0. If there were such 
a root k = e'*,£ real, then the periodic problem would have growing solutions with time. 
This contradicts the fact that the symbol of the operator P~ l Q has purely imaginary 
eigenvalues. The roots k are continuous functions of 3 except at 3 = 3. In the limit, as J 
tends to oo, 

Ki = -2 + \/3 ,|«i| < 1, 

/c 2 = — 2 — \/3, |/c 2 | > 1 • 

Furthermore, «i is continuous at the exceptional point 3 = 3, and we have /Ci = — 1 /2, /c 2 = 
oo at this point. For all other 3 in the right half-plane we have |« 2 | > 1. 

Therefore, if Re(s) > 0, there is only one root k x of (15) with |«i| < 1. and the general 
solution of (12) is 

Uj — 0 \k\. 

Inserting the solution into the boundary equation (13) gives 


D{3)a x = hg , 



where 

D(3) = s(l + 2/cj) + -(5 - 4/cj - /cj). 
The Kreiss condition is fulfilled if 


D(s) 0, Re(s) > 0. 

Assume that this condition is not satisfied. Then we solve the equation D(3) = 0, and 
substitute 

3 = -(—5 -I- 4«i + «i)/(l + 2 ki), K\ ^ — 1/2 

into (15). The resulting equation has the only possible solution /ci = 1, corresponding 
to 3 = 0 in (15). But a perturbation calculation with 3 > 0, s << 1 in (15) shows that 
Kj = — 1 ,k 2 = 1. Since D(3) ^ 0 also for the exceptional value /cj = —1/2, we have shown 
that the Kreiss condition is fulfilled. Therefore we get from (16) 

|<t, | < const. |h< 7 |, 


i.e. 

|uj| < const.|/i£|, Re(s)>0, j = 0,1,.... 

By integrating |ujj 2 along the line Re(s) = 0 and using Parseval’s relation, wo get 


f \uAt)\ 2 dr < const, f \hg(r)\ 2 dr , j ss 0,1,.... 
Jo Jo 
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But Uj ( T ),r < U does not depend on $(r),r > t. Therefore, when considering uAr 
[0, t], we can as well set o(r) — 0 in it. noi Thic _.. 8 J 


;!I set g(r) = 0 in (t, oo). This gives the estimate 
l Mr)| J dr < const, f |hy(r)|‘<fr , .7 = 0 , 1 .... . 

JO 


) > n 


(17) 


Jn/u"wol'T te ‘ Im tai "j d !‘ y T g the energy method - R « alli "* ‘he definition (9) 
O’iQ J)’ e 9- (11)> nnd that P 1 $ SPD. we have for some constants o 


Pu) = 2(u,Qu) + 2uf)hg = T a.,u iU , + 2u„kg , 


*,J=0 


implying 


2 

IMOII 2 < const.(u(t),Pu(t )) < const.{ f (£ h(r)| 2 + \hg{r)\ 2 )dT . 

•'O .• rt 


The final estimate now follows by using (17). 

auxLyml, bZ Pr ° Ve th<! lpi>r0xima,i0 " is stro "*‘y sta hle- Consider first the 


( p -J7h = (Qvh, j = 1 , 2 ,..., 


vo = v,, 


v ;(0) = fj. 


( 18 ) 


Kprun'' 18 th n b .° U n nd Y y co " diti ' ,n wilh r “P«' we can eliminate it */</( such 
t>/ ), is well defined also for J = 1. We now use the scalar product and norm 


(n, t>)i iO0 — ujVjk , 


U l 1,00 = («.«)l,oo. 


and by applying the energy method we obtain 


After integration we get 


dt^ V ' Pv ^ l '°° ~ P v t)t.oo = 2(n, (?r)i t< 

- ~ l ’i u o = -M 2 = -h| 2 . 


h*). Pv(t )) t>00 = (v(0), ^(O)),.* - X - J\\v 0 (t)\ 2 + |t; 1 (r)| 2 )dr . 
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By using the boundary condition, Vo can be eliminated, and it is easily shown that P is 
positive definite. Since P is bounded, we get the estimate 

IM0lli,oo+ / (M r )| 2 + Mr)| 2 )dr < const.\\f\\* t00 . (19) 

j o 

Since the original boundary condition employs 3 points tzo, ui,u 2 , we also need an 
estimate for v 2 . For j = 2,3,... we write (18) as 

dv 

{P Tt )] = {Qv) " > = 2 ' 3 --- 


V\ = Vi, 

«i( 0 ) = 

where uj in the right hand side of the boundary condition is considered as a known 
function. For this new problem we use the same technique as above. We construct a new 
auxiliary problem for which we can derive an energy estimate including the boundary 
values, and for the remaining part of the solution we use the Laplace transform technique 
and the Kreiss condition to obtain an estimate. This time the auxiliary problem is 

._ dw , , ^ 


wi = w 2 , 


“i(0) = />• 

In the same way as we obtained the estimate (19), we now get 

IMOilloo + / (K(t)| 2 + |u> 2 (r)| 2 )dr < cons<.||/||| iOC . (20) 

0 

The grid function y } = v 3 - w 3 satisfies 

( p ^h = (Qy)>, > = 2,3,..., 


!/i = 9 \ i jfi = ui - u>i , 


( 21 ) 


y>(0) = 0 


We have 


Lemma 3.3 The solution of (21) satisfies the estimate 

/ |yj(T)| J dr < const. / |^,(r)| 2 dr, 

Jo 
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i = 1.2 






▼ 


▼ 




Proof: The Laplace transform of (21) is 


~ 2 & +I “&'-»)* s-.sh, j = 2,3,, 


Vi = 9 \, 


llylli.oo < oo 


with the solution 


Vi = 1 , 2 ,... , 


where*, ls the solution of the characteristic equation (15) with |*,| < | for Re(i) > 0 As 

fiTT , ?‘ h r Pr0 : f0f Lemma3 2 . «. U well defined for all i, ke(i) > 0. By il.ratine 
lyj along the line Res = 0, and using ParsevaPs relation, we get ° g 

l Mr)\ 2 dr< const.£°\ gi (r^dr, > = 1,2,. 

follows 16 ”" 113 2 We C “ Cha ” ge l ° inte * ra,i ° n 0VCr a finite «">* interval, and the lemma 

□ 

By this lemma, the definition of o, and the estimatpc 10 m , 

estimate y estimates (iy), (20), we now have an 

l b 2 (r)| 2 dr < const. /‘(|j/ 2 (r)| 2 + |u, 2 (r)| 2 )d T 

/< 

< const. J o (|o,(r)|> + |«,,(r)| ! + |^( T )|t )rfr < const.||/||;^ . 

We also need estimates for dv,/dt near the boundary. By differentiation 1181 with 
res_pect to I we get the same diffe.ential-difference equation and boundary conkll for 

t d ldL t Sl " CC a any t,me *’ we can so,ve boundedly f or dv/dt in terms of Qv we also 
have an initial condition for <f >, yielding the problem ’ 

^ P ~dt ^ = » i~ 2,3,..., 

<Po = <^i, 

m) = \(rs),. 

Here R is a bounded operator, and accordingly, 


1,00 < const.h 


-l 


1,00 


We now use the same procedure as above to derive estimates for A The 
summarized in ^ 


results are 
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Lemma 3.4 The solution of (18) satisfies the estimate 


d-v(t)„ 3 


dt' 


IIU + 


/.5 


d l ‘v } ( T ) l2 


dt 


\dr < const.h 


-2i/ 


1,00 1 V = 0 , 1 . 


□ 

We can now derive the final estimate for u. The difference z 3 = u 3 —v } (where u and 
v are the solutions of (11) and (18) respectively) satisfies 

(p d f t ), = m,, j = 

(P^)o = «Jz)o + 9-(/ > f)o4(^)o, 

-?,(()) = 0. 

The operator P is bounded in the maximum norm, and Q is of the order h~ l . By applying 
Lemma 3.2 and Lemma 3.4, we get 

ll*(OII 2 +/ M t )I 2 dr < const, j (|M T )I 2 + 5Z M T )| 2 + £ \h—jp-\ 2 )dr 

< const .{||/|| 2 + f |^(r)| 2 dr), j = 0.1,.... 

Jo 

By the definition of z and by using Lemma 3.4 once again, we have proved 
Theorem 3.1 The approximation (11) is strongly stable, and the solution satisfies 

IMOII 2 + / £K( t )| 2< * t < consf.(||/|| 2 + [ \hg(T)\ 7 dT). 

Jo Js0 Jo 


□ 

If a forcing function Fj(t) is introduced into the first equation of (11), an extra term 
So ||F(r)|| 2 dr enters the right hand side, see (5). The error estimate then follows immedi¬ 
ately from strong stability. The error e } (t) = u(x,,<) — ufit) satisfies 

(Pj t ), = (Q e )j + 0(h*) , j — 1,2. 

(Pj t ) o = (Q e )o + 0(h 3 ), 

e,(0) = 0, 

and we get from Theorem 3.1 
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Corollary 3.1 The solution of (11) satisfies the error estimate 


|u(xj, t) — Uj(<)|| < const.h 4 . 


□ 


At this point we have finished the analysis of the outflow problem, and we turn to the 
inflow problem 


u t = —u x , 0 < x < oo, 0 < f, 

u(0,f) = g(t), (22) 

u(x,0) = /(*). 


When using the implicit difference operator to compute an approximation Vj of u x {xj, t), 
we need a boundary condition for v r From the differential equation we get, after differ¬ 
entiating the boundary condition with respect to t , 


u r (0,<) = -g'{t), 


and the approximation becomes 


-(Vj-! + 4Vj + v )+i ) 



Uj-i), j = 1,2,..., 


u 0 (t) = flr(f), 
«o(0 = -g'(t), 


which yields 









«o = g, 


(24) 


«,( 0 ) = fj. 

where P is defined at inner points as in (9). Here it is tacitly understood that the 
differentiated form of the boundary condition is used to define du 0 /dt. Clearly, P is SPD 
in the space of grid functions { u ,}" with compact support. 

Corresponding to Lemma 3.2 we have 
Lemma 3.5 Consider (24) with f — 0. The solution satisfies the estimate 
IMOIIi.oo + / K( r )l a< * T < const. I |<7(r)| 2 dr) , j = 1,2,... . 

JO Jo 
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Proof: The Laplace transformed problem is 

s(Pu)i = -(<?«),, j = 1,2,..., 

«o = 9, 

with the characteristic equation 

s(k 2 + 4tc + 1) = ~3(k 2 - 1), s-sh. (25) 

The coefficient of k 2 vanishes for s = -3, which does not cause any trouble, since we are 
only interested in s located in the right half-plane. Therefore, we get immediately 

ii } - gti \, 

where kj is the solution of the characteristic equation (25) with |/ci| < 1 for Re(J) > 0. 
Parseval’s relation yields 

/ \uj(r)\ 2 dr < const, f \g(r)\ 2 dr , j= 0,1,.... (26) 

Jo Jo 

Applying the energy method and using (26) together with the fact that P is SPD proves 
the lemma. 

Remark: For the outflow problem there is a gain of one power of h in the estimate with 
respect to the boundary data. For the inflow problem with a physical boundary condition, 
this gain does not occur. 

In order to prove strong stability, we use the same procedure as for the outflow problem; 
it now becomes much simpler. As our auxiliary problem we now take 

(P—)j = ~(Q v )jy J = l,2, 

Vo = — V\ , 

”;(0) = fj i 

leading to the estimate 

||v(0llico+ / l^o( T )| 2 dr < consf.il/lli oo. (27) 

Jo 

The difference w, = — Vj satisfies 

(p~), = «<»),, j= 1,2 . 

u>o = 9 - Vo, 

u>j(0) = 0. 
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Lemma 3.5 now gives 


IM0lli,oo < con st- f (Mr)| 2 + vo(r)| 2 )dr 

Jo 

< cons<.(||/||| t00 + / \g(T)\ 2 dr ). 

J 0 

The stability follows by using the definition of w and (27): 

Theorem 3.2 The approximation (24) is strongly stable, and the solution satisfies 

IMOIImo ^ C0ns< -(ll/llj.oo + / l$( r )l 2 «fr)- 

Jo 

□ 

The only truncation error occurs in the difference approximation at inner points, and 
an 0(h 4 ) error estimate follows immediately. 

Remark: The method of deriving stability and error estimates presented here can be 
generalized to systems of PDE. For the simple model example treated here we could have 
used the following direct method. □ 

Let 4>{x,t) be a smooth function with 

<f>(x,0) = /(*), 

<£( 0,0 = </( 0 , 

such that 

/■< d l/+t, 6i- t) ft 

Jo ^~dPdt^^ dT ~ con5< (ll/ll + j[ \9{T)\dr), 0 < 1 / + /* < 1 . 

The difference t>j(0 = u,{t) — <j>(x } ,t) satisfies 

(P^h = (Qo)i + p ,, j - 1,2,... , 

v 0 = 0, 
t’>(0) = o, 

where 

/ ||F(r)||,, 00 <fr < cons/.(ll/H,,oo + / |</(r)|dr). 

Jo Jo 

The energy method gives 

2||P ,/ ^ll,.<.jl|P' /1 Hli.oo = ^\\P'>'v\\\„ = j t (v.Pv), m = 2 (v.Pv,), M 

= 2(v, F)i,oo < cons<.||P 1/3 v||i. 00 ||P 1/3 f’||, i00 . 
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w 


After dividing by ||P 1 ^ 2 u||i,oo and integrating, we get 

\\v{t)\U,oo < const. 11111 ,oo < const, f ||/ >1/2 F(r)|| 1 , 0O dr 

V 0 

< const. [ \\F(T)\\ Xi00 dT < const.{\\f\\ x , 00 + I | 0 (r)|dr), 
Jo JO 

which gives us the estimate for u. 


4 Numerical Experiments 

We will now investigate how the previously analyzed difference methods behave in prac¬ 
tice when applied to two different model problems. We have chosen the linear advection 
equation, which will serve as a simple model for contact discontinuities in fluid dynam¬ 
ics, and Burgers’ equation, which is used to study how the schemes treat shocks and 
rarefaction waves. 

Consider the scalar conservation law 


u e + F x - 0, -1 < x < 1 , t > 0 , 

u(x, 0 ) = f(x). 


At the boundary we prescribe u = g if the characteristic is ingoing. We will consider 
different implementations of the flux derivative t\. 

F X =F X , (c-form), . 

F x = (F - G) x -I- G'u x , (e-form), G = - / F(v)dv . (29) 

F x = F'u x , (p-form), M Jo 


The first expression of eq. (29) is the usual conservative form; the second form corresponds 
to the flux vector splitting that results in an energy estimate. The third variant, finally, is 
the primitive form. These forms will lead to numerical methods with different properties. 
It is possible to give a unified presentation by writing 


where 


The flux F = 


F x = ( qF + 0G) X + (7 F' + 6G')u x , 

(30) 

O 

II 

II 

S 

II 

O 

li 

o 

(c-form), 


0=1 0 = — \ 

It 

o 

II 

r- 

(e-form), 

(31) 

o=0 0 = 0 

o 

II 

il 

r~ 

(p-form). 


F(u) is defined by 




F{u) = u 

_ -.7 lr% 

(advection equation), 

/ n_ t ^ i 1mm \ 

(32) 



Thus, the initial-boundary value problem is defined by eqs. (28), (30), (31), (32). 

Next we formulate the semi-discrete problem 

%+(TWar + 0G) + ( 1 r + tCr)D,)) i - W -T)il) i . >= 0,1 .AT, (33) 

Where u = (u 0 . u*) T is the grid function; F = (F 0 ... F N ) T , G = (G 0 . .. G N ) T , where 
i (Uj ) is t e analytic flux evaluated at u : ; the Gj's are defined analogously. The 
operator F is defined by (F’v), = = 0,...,7V (F'( Uj ) is the Jacobian of the 

analytic flux evaluated at u,) with a similar definition of G'. We shall write F'( U j) = F', 
G'iuj) = G'j for brevity. The operator T represents the analytic boundary conditions 
and g contains the boundary data [12], Finally, D is a difference operator approximating 
d/dx to some order of accuracy; D can be either explicit or implicit. Symbolically we 

write D - P Q, where P and Q are local operators. In the case of explicit operators we 
have P = /, ar.d thus D-Q. 

For explicit difference operators the boundary operator T is defined by 


SoUo , j 
(Tu)j — < Uj , j 
, 6\Un , j 


0 , 

1,2,...,7V 
W, 


- 1 , 


where 


Sn= f 0 if F<;>0, _ f 

| 1 otherwise, 1 | 


0 if Fjv < 0, 
1 otherwise * 


The data g is given by g = (gW(t) x...x g^(t)) T , where g <°) and are the analytic 
boundary data. It follows from eq. (33) that the boundary conditions are fulfilled for all 
time if the initial data satisfies the boundary conditions. It is assumed that the type of 
boundary condition remains the same for all time at a given boundary point. One can 
a.ways restart the process, should there be a change at a time t 0 . We point out that the 
c-, e-, and p-forms lead to identical semi-discrete schemes for the advection equation. 

The implicit scheme is formally obtained by setting T = / and by enforcing the 
boundary conditions explicitly. Hence, 


fa] 

dt 


+ (D(aF + SG) + (yF' + 6G')Du) l = 0, ; = 0,1,..., N, 


(34) 


subject to the constraints 


w 0 = p (0) if > 0, u N = gl" if F' n <0. 

In the following we shall confine ourselves to the fourth-order accurate operator D = P~ X Q 
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discussed in sec. 3. For outgoing characteristics at the boundaries we then have 


tto + 2tii, 


j = 0, 


and 


{Pu)j = < -(uj_i +4uj + Uj+i), j = 1,2 ,...,N - 1, 


2u/v_i + u/v * 


— (—5uo + 4ui + U 2 ) i 
2ft 


J =/v, 


j = 0, 


(Qu)_, — < 2^( d" u j+i) ’ 


j = 1,2,... ,7V - 1, 


I ^ (—un-2 — 1 + 5t is) > j — N ■ 

Since the characteristics are assumed to be outgoing (corresponding to the linear outflow 
case, cf. sec. 3), no analytic boundary conditions need to be enforced. Consequently, the 
semi-discrete scheme (34) becomes 


+ tvj + ((ifF 1 + &G')v)j = 0, j = 0,1,...,/V, 


where v and w satisfy 


(Pv) } = (Qu) } , (Pw)j = (Q(aF + 0G))j , j = 0,1,..., N. 

On the other hand, if there are ingoing characteristics at both boundaries, then no bound¬ 
ary modified stencils are needed since one can use the analytic boundary conditions to 
close P and Q. Consider 

(Pv)t = (Qu)j, 


(Pu>), = (Q(aF + 0G)),, 


which is equivalent to 


1 #J ^ 1 1 1 
g(4v,+u 2 ) = ^u 2 -^uo--r 0 , 

^( 4 u>i -I- ujj) = -r(arfj + “ 2 h^°^° • 


(35) 


Suppose that Fq > 0. Then we have the boundary condition u 0 = <7* 0) , which implies 
F 0 = F(g {0) ) and Go = G(g {0) ). Furthermore, v 0 is an approximation of u x (-\,t). Using 


f; gr 

u * - pi F> F'(gW) 
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we obtain 


t>o = — 


9t 


(0) 


F'(gW) • 

Note that this expression is well-defined since F‘(g {0) ) > 0. Similarly, w 0 approximates 

<*F X (- 1, t) + 0G t (- 1, t) = (of' 4 0G')u x (- 1,0 . 


Using 


u, = — 


9i 


( 0 ) 


F'(gM) 

leads to the following boundary approximation for wq 


w 0 


(, . a cV 0> ) \ to 

v +<, FV 0| )J 9 ' 


Substituting these expressions into eq. (35) yields 

— (4l>l 4 O 2 ) = —-U 2 - g^ -\ -—- 

6 V } 2h 2h y 6F'{gl°))’ 


-(4u>, 4 w 2 ) - + 0G 2 ) - —(aF(g {0) ) 4 0G(g {o) )) 4 4 ^ • 

The right boundary is treated analogously. Summing up, the fourth order implicit scheme 
is defined by eqs. (36) - (41) 

du , 

+v>i + {{lF +6G)v)j = 0, j = 0,1,..., JV, (36) 

where v and w are the solutions of 

(Pv) } = (Qu) } + P] , (Pw) } = (Q(aF + 0G))j+q : , j = 0,1,..., N . (37) 

The explicit structure of the difference operators P and Q will be given shortly. We 
have moved the boundary such that x_, = -1 and x/v+i = 1 in eqs. (36) - (41) in 
case of ingoing characteristics. These extra boundary points have then eliminated by 
means of the analytic boundary conditions. This procedure will simplify the computer 
implementation of the algorithm, since the number of unknowns will be same regardless 
of the direction the characteiistics (u_i and u/v +1 are known if the characteristics enter 
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the domain). 
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and 


2h (oF(f< 01 ) + W 1 )) + | (o + 9,"" if Fi 


> 0 , 


0 


j = o. 


otherwise 


<b={ o, 


j = 1,2 ,...,N - 


^ aF(s,,,) + + | (“ ■+ $$$ if Fi < 0, 


0 


j — N, 


otherwise. 


(41) 


The numerical method (36) - (41) i s disrr*ti™t • 

TVD Runge-Kutta method 11.1 ' d llm< ' “ sm * an ex P lici « 3rd order 


[15] 

“ ( " = u"-*£(„"), 

“ l21 = 


(42) 


U 


n+1 _ 


ju n + |ti< 2 >- ~L(u™), 


— < 3 «) 

TVD Runge-Kutta method because of its sim r ^ 'i? S * USC th<? above th «rd-order 

Runge-Kutta methods of higher order of accurJcy'tlfan th’* t0 COnstruct TVD 

more complicated. y ban tbree ’ ^ut they are considerably 

higRorder'centered^ difference a^moxi mat ions T"*?- »hen constructing 

shocks and contact d : scontinuitif»« c ... ,S assum P tlon ' s obviously violated at 

overcome this Z °“ ^ 

:• ,42) As * preii — 

-r = ? + «;+;) = + 1 a+ a_ s -' , > = 1>2 ,. N _, 

a switch r, to ^ ” ln,r ° duC, ‘ 


u; +l = u n+1 


' +^ + (r J . 1/J Aj;+* )t >=1,2. N- 1. 


(43) 


21 


We have used a switch proposed by Jameson [6] 


/ - A-u,| \ m 

\|A + iij| + jA-UjI/ 



1 , 2 - 1 , 



where we have omitted the time index n + 1 in the right member for simplicity. If A + iij 
and A.Uj do not have the same sign, which happens for high-frequent oscillations, then 
r } — 1. For the remaining grid points we obviously have 0 < rj < 1. Taking m = oo 
yields rj = 0 away from the spurious regions. The complete numerical algorithm is thus 
defined by eqs. (36) - (44). 

In the first numerical experiment we solve the linear advection equation to see how well 
the explicit (E) and implicit (I) fourth order methods capture an oscillating solution with 
a discontinuity in the derivative; for second order methods one can expect poor resolution 
at the point of discontinuity. The results are compared with a those of a standard explicit 
second order centered finite difference method and a third order accurate TVD method 
using the following limiter function [15] 


0 



V 


2r 


(2r + 1 )/3 
2 


r < 0, 

0 < r < 1/4 , 
1/4 < r < 5/2, 
r > 5/2. 


As initial data we use 


R(x, a, k) 


( 


—a sin(t7rx) 
x 


x < 0, 
x > 0, 


with a = 0.1 and k = 6. The solutions are plotted at t = 0.5. 


u(x.O) -*<«.0 t.i). u(-1.Q«0. t- 0 5. n.too, CF1.0I6 
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Fig 2. 4th order (E) (+) versus true ( —) solution 



Fig 3. 4th order (I) (-f) versus true ( —) solution 


up.O) ■ Afr.O 1 .0). u(-1. 1 ) • 0, t - 05. n • too. CR • 0 OS 



Fig 4. 3rd order TVD (+) versus true (—) solution 

The fourth order methods clearly resolve the discontinuities much better. It should be 
noted that no artificial viscosity (filter) has been used for the 2nd and 4th order schemes. 














Next we study how many grid points are needed to achieve a certain tolerance level 
£ = II^a — u||<xm where u * is the numerically computed solution. Again we use R(x,0.1,6) 
as initial data and present the solutions at t = 0.5. We have chosen e « 0.04. 


ufr.0)«n(x.01.«). uf0.t)-0, t« OS. n* 100. CFl -006 


Fig 5. 2nd order, ||«a — u||oo = 0.042 for 100 grid points 


u(x.0)»ftQ'.01.t). u(0.t) >0. t.0 5. n-44. CFL-006 


Fig 6. 4th order (E), — uHo, = 0.032 for 48 grid points 


u(k.0)-R{x.61.0) ufO.n-O. t.0(. n.M. CR..006 


0 0 


Fig 7. 4th order (1), ||u^ - u||oo = 0.031 for 34 grid points 


















T 


Thus, the fourth order methods achieve the same level of accuracy as the second order 
method using only half (explicit) or one third (implicit) of the number of grid points. No 
artificial viscosity was used. We have set the CFL-number to 0.05 to suppress errors due 
to the time discretization. 

To simulate contact discontinuities we again solve the linear advection equation, this 
time using piecewise continuous initialdata 


H(x,u L ,u fl ) 


ul x < 0, 
ur x > 0. 


with u L = 1, ur = 0. At x = -1 we prescribe u(-l,<) = 0. The resulting solution is a 
square wave traveling with speed 1 to the right. The fourth order methods are compared 
with the standard second order method and a third order TVD scheme. It is evident 
from the following figures that the fourth order methods are superior to the second order 
method. In fact, the fourth order solutions are comparable to that of the third order TVD 
scheme. For the centered difference schemes we used the previously described filter. 


l 0 

Q 

0 

I. 
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Fig 8. Contact discontinuity, 2nd order 
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Fig 9. Contact discontinuity, 4th order (E) 
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u<xO).H<x1,0), u<-1.t)-0. t.Ot. n-100. Ca*ON 



1 -04 -0* -04 -02 0 02 04 00 01 1 


Fig 10. Contact discontinuity, 4th order (I) 

ufr,0)«H(x.1.0). uM.Q.O. t• 0.5. n.100. CR-0 05 



•1 -00 -Of 4)4 -02 0 02 04 00 00 1 

Fig 11. Contact discontinuity, 3rd order TVD 

We next solve the Riemann problem for Burgers’ equation. For shocks we have used 
the initial data H(x,l,0) and the boundary data u( — l,f) = 1. We include the 3rd order 
TVD solution as a reference. 

up.O)»M0M.O). 1, t.OS n.100 CFl«0tt 



1 *00 40 -44 -01 0 02 04 00 00 I 

Fig 12. Shock, 4th order (E) 
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u(K.O).H0t.1.O). u(*1,t)« 1. t-os. n* 100. CfL.OM 
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♦0» -06 -0.4 -02 6 0.2 04 06 06 1 

Fig 13. Shock, 4th order (I) 
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Fig 14. Shock, 3rd order TVD 


The e-form of the fourth order methods was used, since it appears to be less oscillatory at 
the shock than the other forms. The filter was turned on in a neighborhood of the shock. 
The fourth order methods generate almost as crisp a shock profile - albeit somewhat more 
oscillatory - as the 3rd order TVD scheme. 


Finally we solve Burgers equation for a rarefaction wave. As initial data we take 
H(x, — 1,1). For the fourth order methods we use the c-form as well as the e-form without 
artificial viscosity. The c-form evidently violates the entropy conditions, whereas the 
e-form produces an entropy satisfying rarefaction wave. 
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« * 1 1_A k 1- *- - * i ■ i " - 

• f 4| -01-0402 0 02 04 06 01 1 

Fig 18. Rarefaction wave, 4th order (I), e-form 


U(R.0)«H(X,-1.1) tmOi n• 100. CFL• 066 



Fig 19. Rarefaction wave, 3rd order TVD 


5 Discussion 

In this paper we have studied explicit and implicit fourth order difference operators for 
hyperbolic initial-boundary value problems. Presently there exists no general procedure 
for establishing stability of implicit high-order difference approximations if one wants to 
enforce the analytic boundary conditions explicitly [5], cf. eq. (23). We have presented 
a complete stability analysis of an implicit fourth order accurate difference operator for 
the initial-boundary value problem. The boundary points are eliminated by means of 
the analytic boundary conditions; at boundary points where there is no analytic bound¬ 
ary conditions we use one-sided stencils. The stability result then follows using Laplace 
transform techniques. The result of the stability analysis forms the basis for the actual 
computer implementation of the fourth order implicit operator. 

For implicit and explicit difference operators having one-sided differences at every 
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boundary point there is a well-developed stability theory based on the energy method. 
The analytic boundary conditions are then enforced by adding a penalty function or a 
projection to the semi-discrete system [3, 12]. We have followed the approach in [12] for 
the implementation of the explicit fourth order operator . 

It has been numerically verified that the fourth order methods studied in this paper are 
more efficient than the standard second order one. For the linear test problem, figures 1 - 
4 show that discontinuities in the derivative and high frequency data are better resolved 
using the high-order difference methods. Also, they are more efficient to achieve a certain 
tolerance level (figures 5-7). In the one-dimensional case we obtained a reduction of grid 
points by a factor of two for the explicit fourth order method, and by a factor of three for 
the implicit operator. This is true for each space dimension. Thus, in three dimensions 
one would obtain a reduction by a factor of eight or twenty-seven, respectively. Since 
the work grows linearly it is natural to assume that high-order methods would be even 
more efficient for multidimensional problems. We emphasize that no artificial viscosity 
was used in the previous test cases. 

The fourth order methods are good candidates for handling the case where the data is 
piecewise continuous. This is illustrated in figures 8 - 11 Artificial viscosity was needed 
to control spurious oscillations in this case. The performance of the fourth order methods 
is comparable to that of a third order TVD method. 

The numerical experiments were concluded by solving Burgers’ equation. Two different 
forms of the flux derivative was implemented: the c-form and the e-form. The c-form is 
the usual conservative form, and it may lead to entropy violating solutions for both the 
implicit and explicit operators, see figures 15 and 17. The e-form, however, picked up the 
entropy satisfying solution without using artificial viscosity (figures 16 and 18). Indeed, 
in a forthcoming paper it will be shown that for diagonal norms H one can prove an 
entropy condition for the semi-discrete system if the e-form is used. Furthermore, shocks 
are treated satisfactorily after adding artificial viscosity, cf. figures 12 - 14. 

In summary, there is a complete stability theory for the high-order methods that 
we have used. The theoretical properties have been verified through numerical experi¬ 
ments. For nonlinear conservation laws these high-order methods work as well as specially 
constructed high-order TVD schemes; for linear problems with high frequency solutions 
(or discontinuities in the derivative) the difference methods work better than the TVD 
schemes. Another attractive feature of these difference methods is the simplicity of their 
computer implementation. We anticipate that these methods will generalize well to sys¬ 
tems of conservation laws, where all phenomena (shocks, contact discontinuities, etc.) 
may be present at the same time. 
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6 Appendix 


Fourth order accurate difference operator with third order boundary closure satisfying 
eq. (5): 


{ 1 


— (dooUo + doiUi + do2U2 + do 3 « 3 ) 

h 


j = 0 


(Du) } = 


Y(dioUo ■+■ duuj + di2U2 + di3«3 + dj4«4 + di 5 Us) j — 1 

/i 

-)-(d20Uo + djjUi + d22U2 + ^23^3 + ^ 24^4 + <^ 25 ^ 5 ) J = 2 

h 

t(^ 30«0 + d3iUj + dyiV ,2 + d 33«3 + <^34114 + + d^U^) j= 3 

n 

y(c< 4 oUo + (/ 41 U 1 + d 4 2U2 + d 43 ti3 + d 44 U4 + d 45 u 5 + d 4 eu e ) j = 4 

h 


\2h 


(uj_2 — 8 Uj_i + 8uj+i — Uj+ 2 ) 


j = 5,6, 


The corresponding norm is defined by 


(Hu) } = 


( hoouo j 

h\\U\ + h\2U2 + h\3U3 + /»14 u 4 j 
h\2Ui + /l22 u 2 + /l23 l *3 "1" ^24^4 j 
' h\3U\ + /l23U2 + A 33 U 3 + ^34 u 4 j 
J*t4«l + ^24“2 + ^ 34^3 + ^44^4 j 


0 

1 

2 

3 

4 




j — 5,6, • • • 


The elements dij are given by 


doo 

doi 

do2 

do3 



f\d\o 

Sxdn 

f\du 

fidi3 
/id ,4 
/idis 


-24(-779042810827742869 + 104535124033147\/26116897) 
-(-176530817412806109689 + 297682748168?5927 >/2611689 7) /6 
343(-171079116122226871 + 27975630462649^ 26116897) 

- 3(- 7475554291248533227 + 1648464218793925\/26116897 )/2 
(-2383792768180030915 + 1179620587812973> /26116897 )/3 
— 1232(—115724529581315 + 37280576429\/26116897) 
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fid 20 
fid 2 \ 
fid 22 
fld 23 
fid 24 
fld 2 s 


— 12(—380966843 + 86315/26116897) 
(5024933015 + 2010631/26116897 )/3 
—231(—431968921 + 86711/26116897 )/2 
(-65931742559 + 12256337/26116897) 
-(-50597298167 + 9716873/26116897 )/6 
-88(-15453061 + 2911/26116897) 


fxdao = 48(—56020909845192541 + 9790180507043\/26116897) 

fid 3 ! = (-9918249049237586011 + 1463702013196501/26116897 )/6 

fid 32 = —13(—4130451756851441723 + 664278707201077/26116897) 

/W 33 = 3(-26937108467782666617 + 5169063172799767/26116897 )/2 
fxd 34 = -(6548308508012371315 + 3968886380989379/26116897 )/3 

f\d 3 $ = 88(—91337851897923397 + 19696768305507/26116897) 

Mk, = 242(-120683 + 15/26116897) 



264(-120683 + 15/26116897) 

(-43118111 + 23357/26116897 )/3 
—47(-28770085 + 2259/26116897 )/2 
-3(1003619433 + 11777/26116897) 

—11(—384168269 + 65747/26116897 )/6 
22(87290207 + 10221/26116897) 
-66(3692405 + 419/26116897) 

-56764003702447356523 + 8154993476273221 /26116897 
-55804550303 + 9650225/26116897 
3262210757 + 271861/26116897 


The elements h{j are j'iven by 


hoo = 3/11 

fh u = (299913292801 + 56278767/26116897)/228096 
fh u = -(64756272879 + 310129/2611 €C97)/76032 

fh\ 3 = -(-50615837729 + 5284177/26116897 )/76032 
fh u = (-5026701941 + 948741/26116897)/20736 

fh 22 = -7(-6989673895 + 13527/26116897 )/25344 

jh 23 - 49(—657605303 + 100423/26116897)/25344 

fh 24 = —49(—75022899 + 14467/26116897 )/6912 

//»33 = -(-45333081425 + 982369/26116897 J./25344 

//»34 “ (-3355209517 + 597005/26116897)/6912 

fh 44 = 5(35213725709 + 5139171/26116897 )/228096 
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where/ = 591223+ 146^26116897. In decimal form the elements d„ can be expressed as 

doo = -1 83333333333333333333333333333 

do i = 3 

d 02 = -1 50000000000000000000000000000 

d 0 3 = 0.333333333333333333333333333333 

dio = -0.389422071485311842975177265601 

d u = -0.269537639034869460503559633382 

du = 0.639037937659262938432677856177 

d n = 0.0943327360845463774750968877542 

du = -0.0805183715808445133581024825053 
d xs = 0.00610740835721650092906463755986 

<*20 = 0.111249966676253227197631191910 

<*21 = -0.786153109432785509340645292043 

d 2 2 = 0.198779437635276432052935915731 

<*23 = 0.508080676928351487908752085978 

<*24 = -0.0241370624126563706018867104972 
<*25 = -0.00781990939443926721678719106473 

<*30 = 0.0190512060948850190478223587424 

<*3i = 0.0269311042007326141816664674714 

<*32 = -0.633860292039252305642283500160 

<*33 = 0.0517726709186493664626888177642 

<*34 = <1.592764606048964306931634491846 

<*35 = -0.0543688142698406758774679261364 
<*36 = -0.00229048095413832510406070952285 

<* 4 o = -0.00249870649542362738624804675220 
<*4i = 0.00546392445304455008494236684033 
<*42 = 0.0870248056190193154450416111555 

<*43 = -0.686097670431383548237962511317 

<*44 = 0.0189855304809436619*79348998897 

<*45 = 0.659895344563505072850627735852 

<*46 = -0.0827732281897054247443360556719 
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